Bur Write Up

Data Context & Research Question

We will examine two separate tree growth data sets with the same variables. Both data sets, bur and pin, contain a tree id number, location data, the initial diameter of the tree (measured at breast height), logged soil nitrogen levels, the elevation of the tree, and its growth over 25 years. We will investigate whether or not the elevation of a tree will affect its growth while taking into account different tree diameters. We will also run a similar analysis and see whether or not the logged soil nitrogen levels affect a tree’s growth while taking each tree’s starting diameter. Both analysis’ will be run for the bur and pin dataset resulting in four total models.

Visualizations

Tree Maps

bur %>%
    ggplot(aes(x = X, y = Y, color = `ELEVA`)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Bur Trees Map with Elevation", color = "Elevation")

bur %>%
    ggplot(aes(x = X, y = Y, color = `LOG(N)`)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Bur Trees Map with Logged Soil Nitrogen Levels", color = "Logged Soil Nitrogen")

pin_growth %>%
    ggplot(aes(x = X, y = Y, color = `ELEVA`)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Pin Trees Map with Elevation", color = "Elevation")

pin_growth %>%
    ggplot(aes(x = X, y = Y, color = `LOG(N)`)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Pin Trees Map with Logged Soil Nitrogen Levels", color = "Logged Soil Nitrogen")

trees %>%
    ggplot(aes(x = X, y = Y, color = `ELEVA`,shape = type)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Trees Map with Elevation", color = "Elevation" ,shape = 'Tree Type')

trees %>%
    ggplot(aes(x = X, y = Y, color = `LOG(N)`,shape = type)) + 
    geom_point() + 
    theme_classic() +
    labs(title = "Trees Map with Logged Soil Nitrogen Levels", color = "Logged Soil Nitrogen",shape = 'Tree Type')

Our initial investigations show that towards the left side of the visualization, where X = 0 to X = 100, elevation is lowest. At X = 200, Y = 200-400, there seems to be a hill as elevation is around ~284. We can also see another hill in the bottom right of the map. The logged soil nitrogen levels are almost the same across the entire map besides a small cluster on the left-hand side at X = 50, Y = 175.

Examining Growth

trees %>%
  ggplot(aes(x = GROWTH, fill = type)) +
  geom_density(alpha=.2) +
  labs(x = "Growth", title = "Distribution of Growth") +
  theme_classic() +
    theme(plot.title = element_text(hjust = 0.5)) 

trees %>% count(type, GROWTH >= 0) # 5 bur trees & 3 pin trees with negative growth

##   type GROWTH >= 0    n
## 1  bur       FALSE    5
## 2  bur        TRUE 2238
## 3  pin       FALSE    3
## 4  pin        TRUE  441

trees %>% filter(GROWTH < 0)

##     ID_     X     Y DBH1995 GROWTH   LOG(N)  ELEVA type
## 1 11790 208.0 189.9    15.0  -15.0 6.753438 282.92  pin
## 2 11794 377.0 184.0    14.5  -14.5 6.705639 283.69  pin
## 3 11838 353.0 241.0    13.8  -13.8 6.448889 282.52  pin
## 4  1147 113.0  87.0    48.0  -48.0 6.526495 281.11  bur
## 5  1542  92.0 108.0     6.3  -14.6 6.639876 280.63  bur
## 6  6181  65.5  41.0     8.5   -8.5 6.876265 280.25  bur
## 7  7113 159.8 177.0     8.5   -8.5 6.870053 282.62  bur
## 8  9359 240.1  61.0     7.0   -7.0 6.834109 280.52  bur

trees %>% filter(GROWTH >= 0) %>%
  ggplot(aes(x = GROWTH, fill = type)) +
  geom_density(alpha = .2) +
  labs(x = "Growth", title = "Distribution of Growth (> 0) among Bur Trees",fill = 'Tree Type') +
  theme_classic() +
    theme(plot.title = element_text(hjust = 0.5)) 

Across all exploratory visualizations, we find five bur and three pin trees with negative growth.

Examining Variable Relationships

trees %>% filter(GROWTH >= 0) %>%
  ggplot(aes(x = DBH1995, y = GROWTH, color = type)) +
  geom_point(alpha = .2) +
  geom_smooth(se=FALSE) + 
  labs(x = "Diameter of Trunk at Breast Height", y = "Growth", title = "Trunk Diameter vs Growth",  color = 'Tree Type') +
  theme_classic() +
   theme(plot.title = element_text(hjust = 0.5))

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

p1 <- trees %>% filter(GROWTH >= 0) %>% 
  ggplot(aes(x = `LOG(N)`, y = GROWTH, color = type)) +
  geom_point(alpha = 0.2) +
  geom_smooth(se=FALSE) + 
  labs(x = "Logged Soil Nitrogen Levels", y = "Growth", title = "Logged Soil Nitrogen vs Growth",  color = 'Tree Type') +
  theme_classic() +
   theme(plot.title = element_text(hjust = 0.5))

p2 <- trees %>% filter(GROWTH >= 0) %>%
  ggplot(aes(x = ELEVA, y = GROWTH, color = type)) +
  geom_point(alpha = 0.2) +
  geom_smooth(se=FALSE) + 
  labs(x = "Elevation", y = "Growth", title = "Elevation vs Growth",color = 'Tree Type') +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) 

grid.arrange(p1, p2, ncol = 2)

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

We can see that the majority of trees have an trunk diameter between 5 and 20 units. Also, when looking at the logged soil nitrogen levels, we see just a few trees have high levels while the majority remain around 5.5 to 7.5, similar to what we saw from the soil nitrogen map above. The elevation of the trees looks relatively normal.

Modeling

Linear Models

trees <- trees %>% filter(GROWTH > 0)

Creating and Interpretating the Models

lm_bur_eleva <- lm(GROWTH ~ ELEVA + DBH1995, data = trees %>% filter(type == 'bur'))

summary(lm_bur_eleva)

## 
## Call:
## lm(formula = GROWTH ~ ELEVA + DBH1995, data = trees %>% filter(type == 
##     "bur"))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.760 -1.693 -0.355  1.274 14.972 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 60.757007   9.341982   6.504 9.64e-11 ***
## ELEVA       -0.203263   0.033191  -6.124 1.08e-09 ***
## DBH1995      0.070684   0.006173  11.451  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.381 on 2235 degrees of freedom
## Multiple R-squared:  0.0679, Adjusted R-squared:  0.06706 
## F-statistic:  81.4 on 2 and 2235 DF,  p-value: < 2.2e-16

Our first model, we will utilize ordinary least-squares to measure how elevation and the logged base diameter of a tree affect the trees growth of bur trees. The intercept does not make sense to interpret in the data context as a tree cannot have a base trunk diameter of zero. We’d expect the average tree to shrink 0.20 units for every 1 unit increase in elevation while holding the base diameter of a tree constant. Additionally, we’d expect the average tree to grow around 0.07 units for every 1 unit increase in base diameter of a tree while holding elevation constant. All the coefficients have statistically significant p-values, but the standard errors maybe invalid due to the spatial correlation of the observations.

lm_pin_eleva <- lm(GROWTH ~ ELEVA + DBH1995, data = trees %>% filter(type == 'pin'))

summary(lm_pin_eleva)

## 
## Call:
## lm(formula = GROWTH ~ ELEVA + DBH1995, data = trees %>% filter(type == 
##     "pin"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.7323  -3.9072  -0.1475   3.2840  17.0021 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 234.76115   48.97370   4.794 2.25e-06 ***
## ELEVA        -0.79101    0.17384  -4.550 6.95e-06 ***
## DBH1995       0.10634    0.02053   5.180 3.39e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.413 on 438 degrees of freedom
## Multiple R-squared:  0.09974,    Adjusted R-squared:  0.09563 
## F-statistic: 24.26 on 2 and 438 DF,  p-value: 1.015e-10

Similarly for pin trees, we will utilize ordinary least-squares to measure how elevation and the logged base diameter of a tree affect the trees growth of bur trees. We’d expect the average tree to shrink 0.80 units for every 1 unit increase in elevation while holding the base diameter of a tree constant. Additionally, we’d expect the average tree to grow around 0.1 units for every 1 unit increase in base diameter of a tree while holding elevation constant. All the coefficients have statistically significant p-values, but the standard errors maybe invalid due to the spatial correlation of the observations.

lm_bur_n <- lm(GROWTH ~ `LOG(N)` + DBH1995, data = trees %>% filter(type == 'bur'))

summary(lm_bur_n)

## 
## Call:
## lm(formula = GROWTH ~ `LOG(N)` + DBH1995, data = trees %>% filter(type == 
##     "bur"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8705 -1.6922 -0.3252  1.3046 14.0478 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.043462   0.700845  -2.916  0.00358 ** 
## `LOG(N)`     0.835276   0.103648   8.059 1.24e-15 ***
## DBH1995      0.070814   0.006134  11.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.366 on 2235 degrees of freedom
## Multiple R-squared:  0.07902,    Adjusted R-squared:  0.0782 
## F-statistic: 95.88 on 2 and 2235 DF,  p-value: < 2.2e-16

In our second ordinary least squares model for bur trees, we measure how logged soil nitrogen levels and the base diameter of the trunk affect tree growth. We’d expect the average tree to grow by 0.831 units for when soil nitrogen levels increase by a factor of 2.7 while holding the base diameter of a tree constant. Additionally, we’d expect the average tree to grow around 0.07 units for every 1 unit increase in base diameter of a tree while holding soil nitrogen levels constant. Also, all the coefficients are statistically significant, but the standard errors maybe invalid due to the spatial correlation of the observations.

lm_pin_n <- lm(GROWTH ~ `LOG(N)` + DBH1995, data = trees %>% filter(type == 'pin'))

summary(lm_pin_n)

## 
## Call:
## lm(formula = GROWTH ~ `LOG(N)` + DBH1995, data = trees %>% filter(type == 
##     "pin"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1344  -3.9278  -0.4058   3.5173  17.6659 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.96420    3.86919  -0.249 0.803323    
## `LOG(N)`     1.90708    0.56637   3.367 0.000826 ***
## DBH1995      0.11338    0.02079   5.453 8.28e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.469 on 438 degrees of freedom
## Multiple R-squared:  0.08097,    Adjusted R-squared:  0.07678 
## F-statistic:  19.3 on 2 and 438 DF,  p-value: 9.305e-09

We’d expect the average pin tree to grow by 1.9 units for when soil nitrogen levels increase by a factor of 2.7 while holding the base diameter of a tree constant. Additionally, we’d expect the average tree to grow around 0.11 units for every 1 unit increase in base diameter of a tree while holding soil nitrogen levels constant. Also, all the coefficients are statistically significant, but the standard errors maybe invalid due to the spatial correlation of the observations.

Limitations of Linear Models

While these linear models are a good starting point, they fail to take into consideration the tree’s location. If two trees are close together, they might be spatially correlated with each other. Without taking this spatial correlation into consideration, we will get inaccurate standard errors for our models, however the slope estimates are still accurate. To fix this problem, we will utilize a spatial mixed-effects model that will account for the location of the trees while also allowing for each tree to have a different intercept for its trunk diameter.

lm.model.map %>% 
  ggplot(aes(x = X, y = Y, color = lm_bur_n_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Soil Nitrogen Model - Bur Trees", color = "Residuals")

lm.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = lm_pin_n_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Soil Nitrogen Model - Pin Trees", color = "Residuals")

lm.model.map %>% 
  ggplot(aes(x = X, y = Y, color = lm_bur_eleva_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Elevation Model - Bur Trees", color = "Residuals")

lm.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = lm_pin_eleva_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Elevation Model - Pin Trees", color = "Residuals")

Spatial Mixed-Effects Models

Quick Overview of spaMM Models

Similar to the ordinary least squares model, for the spatial mixed-effects models I will create two models. One of the models will have elevation as a predictor and the other will have logged soil nitrogen levels. For both the models I will utilize a Matérn covariance matrix. The Matérn covariance matrix is a common method when computing spatial mix-effects models as it will allow us to calculate the Matérn covariance between two points separated by a distance, d. For both models, I will use an unset \(\nu\) parameter to have the optimized values of \(\nu\) and \(\rho\) which produce the lowest residuals. For larger data sets, we could utilize a \(\nu = 0.5\) to make our model more computationally efficient at the cost of accuracy. However, the bur tree data set is small enough where our model with an unset parameter should only take a few minutes to run.

Spatial Elevation Models

spa_bur_eleva_gaussian <- fitme(GROWTH ~ ELEVA + DBH1995 + Matern(1 | X + Y), data = trees %>% filter(type == 'bur'), family = "gaussian")

## (One-time message:) Choosing matrix methods took 5.73 s.
##   If you perform many similarly costly fits, setting the method
##   by control.HLfit=list(algebra=<"spprec"|"spcorr"|"decorr">) may be useful,
##   see help("algebra"). "decorr" has been selected here.

spa_pin_eleva_gaussian <- fitme(GROWTH ~ ELEVA + DBH1995 + Matern(1 | X + Y), data = trees %>% filter(type == 'pin'), family = "gaussian")

spaMM.model.map <- cbind(trees %>% filter(type == 'bur'), spa_bur_eleva_gaussian_pred = predict(spa_bur_eleva_gaussian)) %>% 
    mutate(spa_bur_eleva_gaussian_resid = GROWTH - spa_bur_eleva_gaussian_pred)

spaMM.model.map2 <- cbind(trees %>% filter(type == 'pin'), spa_pin_eleva_gaussian_pred = predict(spa_pin_eleva_gaussian)) %>% 
    mutate(spa_pin_eleva_gaussian_resid = GROWTH - spa_pin_eleva_gaussian_pred)

# spaMM Elevation RMSE
sqrt(mean(spaMM.model.map$spa_bur_eleva_gaussian_resid^2))

## [1] 0.9426431

sqrt(mean(spaMM.model.map2$spa_pin_eleva_gaussian_resid^2))

## [1] 0.4082684

Spatial Soil Nitrogen Models

spa_bur_n_gaussian <- fitme(GROWTH ~ `LOG(N)` + DBH1995 + Matern(1 | X + Y), data = trees %>% filter(type == 'bur'), family = "gaussian")

spa_pin_n_gaussian <- fitme(GROWTH ~ `LOG(N)` + DBH1995 + Matern(1 | X + Y), data = trees %>% filter(type == 'pin'), family = "gaussian")

spaMM.model.map <- spaMM.model.map %>% mutate(spa_bur_n_gaussian_pred = as.numeric(predict(spa_bur_n_gaussian))) %>% 
    mutate(spa_bur_n_gaussian_resid = GROWTH - spa_bur_n_gaussian_pred)

spaMM.model.map2 <- spaMM.model.map2 %>% mutate(spa_pin_n_gaussian_pred = as.numeric(predict(spa_pin_n_gaussian))) %>% 
    mutate(spa_pin_n_gaussian_resid = GROWTH - spa_pin_n_gaussian_pred)

OLS vs Spatial Elevation Model Comparison

Elevation Evaluation Metrics

# OLS RMSE
sqrt(mean(lm.model.map$lm_bur_eleva_resid^2))

## [1] 2.378937

# spaMM RMSE
sqrt(mean(spaMM.model.map$spa_bur_eleva_gaussian_resid^2))

## [1] 0.9426431

We find that our spatially mixed effect model has a significantly lower root mean squared error (0.96) than the OLS elevation model (2.38) for bur trees.

# OLS RMSE
sqrt(mean(lm.model.map2$lm_pin_eleva_resid^2))

## [1] 5.394231

# spaMM RMSE
sqrt(mean(spaMM.model.map2$spa_pin_eleva_gaussian_resid^2))

## [1] 0.4082684

We find that our spatially mixed effect model has a significantly lower root mean squared error (0.40) than the OLS elevation model (5.39) for bur trees.

Elevation Residual Maps

lm.model.map %>% 
  ggplot(aes(x = X, y = Y, color = lm_bur_eleva_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Elevation Model - Bur Trees", color = "Residuals")

spaMM.model.map %>% 
  ggplot(aes(x = X, y = Y, color = spa_bur_eleva_gaussian_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "Spatial Mixed-Effects Elevation Model - Bur Trees", color = "Residuals")

Looking at the residual maps for our OLS and spaMM elevation models we can see that the scale on the right-hand side is much larger for the OLS model. This indicates that the magnitude for the residuals of the spatial mixed-effects model is much smaller, which is ideal. Additionally, when looking at the residual map for the spatial mixed effects model, there does not seem to be any specific area where our model is systematically over or under predicting the residuals. This leads us to believe that there might be no spatial correlation among the tress.

lm.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = lm_pin_eleva_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Elevation Model - Pin Trees", color = "Residuals")

spaMM.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = spa_pin_eleva_gaussian_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "Spatial Mixed-Effects Elevation Model - Pin Trees", color = "Residuals")

OLS vs Spatial Soil Nitrogen Model Comparison

Soil Nitrogen Evaluation Metrics

# OLS RMSE
sqrt(mean(lm.model.map$lm_bur_n_resid^2))

## [1] 2.364703

# spaMM RMSE
sqrt(mean(spaMM.model.map$spa_bur_n_gaussian_resid^2))

## [1] 0.9308835

Similar to our elevation models, the spatially mixed-effects model produce a lower root mean squared error (0.93) compared to the OLS model’s RMSE (2.36) for bur trees.

# OLS RMSE
sqrt(mean(lm.model.map2$lm_pin_n_resid^2))

## [1] 5.450167

# spaMM RMSE
sqrt(mean(spaMM.model.map2$spa_pin_n_gaussian_resid^2))

## [1] 0.4084727

The spatially mixed-effects model produce a lower root mean squared error (0.41) compared to the OLS model’s RMSE score (5.45).

Soil Nitrogen Residual Map

lm.model.map %>% 
  ggplot(aes(x = X, y = Y, color = lm_bur_n_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Soil Nitrogen Model - Bur Trees", color = "Residuals")

spaMM.model.map %>% 
  ggplot(aes(x = X, y = Y, color = spa_bur_n_gaussian_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "Spatial Mixed-Effects Soil Nitrogen Model - Bur Trees", color = "Residuals")

lm.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = lm_pin_n_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "OLS Soil Nitrogen Model - Pin Trees", color = "Residuals")

spaMM.model.map2 %>% 
  ggplot(aes(x = X, y = Y, color = spa_pin_n_gaussian_resid)) +
  geom_point() + 
  theme_classic() +
  scale_color_gradient2(low = "blue", mid = "lightgrey", high = "red") +
  labs(title = "Spatial Mixed-Effects Soil Nitrogen Model - Pin Trees", color = "Residuals")

Again, looking at our residual maps, we see from the scale that the magnitude of the residuals from our spatially mixed effects model are much lower. However, there seems to be no pattern for which areas are over and under predicted, so our data might not be spatially correlated.

Model Summary

Elevation

Using the spatial mixed effects model, there is some evidence that elevation has a negative relationship with growth after accounting for tree diameter. This relationship is strong among pin trees with a -0.67 coefficient that is statistically discernible from 0. Bur trees display a weaker relationship that is moderately statistically discernible from 0.

summary(spa_bur_eleva_gaussian)

## formula: GROWTH ~ ELEVA + DBH1995 + Matern(1 | X + Y)
## ML: Estimation of corrPars, lambda and phi by ML.
##     Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing p_v.
## family: gaussian( link = identity ) 
##  ------------ Fixed effects (beta) ------------
##             Estimate  Cond. SE t-value
## (Intercept) 49.63478 23.427523   2.119
## ELEVA       -0.16410  0.083114  -1.974
## DBH1995      0.09178  0.006154  14.914
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.nu      1.rho 
## 0.12435150 0.02575732 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    X + Y  :  4.008  
## # of obs: 2238; # of groups: X + Y, 2229 
##  -------------- Residual variance  ------------
## phi estimate was 2.01981 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4965.637

summary(spa_pin_eleva_gaussian)

## formula: GROWTH ~ ELEVA + DBH1995 + Matern(1 | X + Y)
## ML: Estimation of corrPars, lambda and phi by ML.
##     Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing p_v.
## family: gaussian( link = identity ) 
##  ------------ Fixed effects (beta) ------------
##             Estimate Cond. SE t-value
## (Intercept) 202.4435 64.72660   3.128
## ELEVA        -0.6761  0.22981  -2.942
## DBH1995       0.1214  0.02128   5.706
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.nu      1.rho 
## 0.05959285 0.03829179 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    X + Y  :  27.22  
## # of obs: 441; # of groups: X + Y, 440 
##  -------------- Residual variance  ------------
## phi estimate was 2.08997 
##  ------------- Likelihood values  -------------
##                       logLik
## p_v(h) (marginal L): -1358.7

Soil Nitrogen

Using the spatial mixed effects model, there is weak evidence that soil nitrogen has a positive relationship with growth after accounting for tree diameter. This relationship is stronger among pin trees with a 1.4 coefficient that is moderately discernible from 0. Bur trees display a weaker relationship that is not statistically discernible from 0.

summary(spa_bur_n_gaussian)

## formula: GROWTH ~ `LOG(N)` + DBH1995 + Matern(1 | X + Y)
## ML: Estimation of corrPars, lambda and phi by ML.
##     Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing p_v.
## family: gaussian( link = identity ) 
##  ------------ Fixed effects (beta) ------------
##             Estimate Cond. SE t-value
## (Intercept)   2.3707 1.220218  1.9428
## `LOG(N)`      0.1527 0.180516  0.8457
## DBH1995       0.0906 0.006131 14.7778
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.nu      1.rho 
## 0.12198483 0.02444717 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    X + Y  :  4.068  
## # of obs: 2238; # of groups: X + Y, 2229 
##  -------------- Residual variance  ------------
## phi estimate was 1.9958 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -4967.19

summary(spa_pin_n_gaussian)

## formula: GROWTH ~ `LOG(N)` + DBH1995 + Matern(1 | X + Y)
## ML: Estimation of corrPars, lambda and phi by ML.
##     Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing p_v.
## family: gaussian( link = identity ) 
##  ------------ Fixed effects (beta) ------------
##             Estimate Cond. SE t-value
## (Intercept)    2.393  4.96857  0.4816
## `LOG(N)`       1.425  0.72977  1.9525
## DBH1995        0.123  0.02147  5.7285
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.nu      1.rho 
## 0.05532534 0.02813536 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    X + Y  :  27.74  
## # of obs: 441; # of groups: X + Y, 440 
##  -------------- Residual variance  ------------
## phi estimate was 2.09623 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1360.478

Conclusion